﻿namespace NCM_MSTest.Alg
{
    /// <summary>
    /// 分段插值
    /// </summary>
    public class PicewiseInterpolation
    {
        public static void LinearInterp(int n, double[] x, double[] y, int m, double[] x1, ref double[] y1)
        {
            for (int m1 = 0; m1 < m; m1++)
            {
                int i = 1;
                for (; i < n; i++)
                {
                    if (x1[m1] <= x[i])
                        break;
                }
                y1[m1] = (x1[m1] - x[i]) / (x[i - 1] - x[i]) * y[i - 1]
                    + (x1[m1] - x[i - 1]) / (x[i] - x[i - 1]) * y[i];
            }
        }

        public static void ThreeInterpA(int n, double[] x, double[] y, int m, double[] x1, ref double[] y1)
        {
            for (int m1 = 0; m1 < m; m1++)
            {
                int i1 = 1;
                for (; i1 < n; i1++)
                {
                    if (x1[m1] <= x[i1])
                        break;
                }
                int i = 0;
                if (i1 != 2 && x[i1] - x1[m1] > x1[m1] - x[i1 - 1])
                {
                    i = i1 - 1;
                }
                else if (i1 == n)
                {
                    i = n - 1;
                }
                else
                {
                    i = i1;
                }
                System.Diagnostics.Trace.WriteLine($"A-->{i}");
                y1[m1] = 0;
                for (int j = i - 1; j <= i + 1; j++)
                {
                    double h = 1;
                    for (int k = i - 1; k <= i + 1; k++)
                    {
                        if (j != k)
                            h = h * (x1[m1] - x[k]) / (x[j] - x[k]);
                    }
                    y1[m1] = y1[m1] + h * y[j];
                }
            }
        }

        public static void ThreeInterpB(int n, double[] x, double[] y, int m, double[] x1, ref double[] y1)
        {
            for (int m1 = 0; m1 < m; m1++)
            {
                int i = n - 2;
                for (int i1 = 1; i1 < n - 2; i1++)
                {
                    if (x1[m1] <= x[i1])
                    {
                        if (i1 != 2
                            && i1 != n - 1
                            && x[i1] - x1[m1] > x1[m1] - x[i1 - 1])
                            i = i1 - 1;
                        else if (i1 == n - 1)
                            i = i1 - 2;
                        else
                            i = i1;
                        break;
                    }
                }
                System.Diagnostics.Trace.WriteLine($"B-->{i}");
                y1[m1] = 0;
                for (int j = i - 1; j <= i + 1; j++)
                {
                    double h = 1;
                    for (int k = i - 1; k <= i + 1; k++)
                    {
                        if (j != k)
                            h = h * (x1[m1] - x[k]) / (x[j] - x[k]);
                    }
                    y1[m1] = y1[m1] + h * y[j];
                }
            }
        }
    }
}
